package fem2.observer;

import inf.common.math.LinAlg;

import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.Set;
import java.util.TreeMap;
import java.util.TreeSet;

import math2.FunctionBasisOnRnToR;
import math2.FunctionRnToR;
import math2.MultFunctionRnToR;
import math2.delaunay.AbstractDelaunay;
import math2.delaunay.Delaunay2;
import math2.delaunay.Triangle;
import fem2.BasisManager;
import fem2.Debugger;
import fem2.Element;
import fem2.IntegrationManager;
import fem2.IntegrationPoint;
import fem2.MeshPart;
import fem2.Model;
import fem2.Node;
import fem2.Observer;
import fem2.enu.EnrichmentType;
import fem2.enu.IntegrationType;
import fem2.level_set.Crack;
import fem2.level_set.Discontinuity;

/**
 * Xfem method for structural problem
 * 
 * @author hbui
 * 
 */
public class XfemManager implements Observer {

	private ArrayList<Crack> crackList;
	/*
	 * TODO: extend for inclusion and void here
	 */

	private Model m;

	private Map<Node, Map<Discontinuity, EnrichmentType>> enrichmentMap = new TreeMap<Node, Map<Discontinuity, EnrichmentType>>();

	int p_on_sub_triangle = 5; // order of integration on sub-triangles

	public XfemManager(Model m) {
		this.m = m;
		crackList = new ArrayList<Crack>();
	}

	/**
	 * add a crack
	 * 
	 * @param c
	 */
	public void addCrack(Crack c) {
		crackList.add(c);
	}

	@Override
	public void fireStarted() {
	}

	@Override
	public void fireTimeStepStarted() {

		/************************************************************************
		 * enrich elements
		 ***********************************************************************/
		/*
		 * build enrichment mappping. Enrichment mapping indicates which type of
		 * enrichment at a node w.r.t each type of discontinuites
		 */
		Set<Element> fullyEnrichedElements = new TreeSet<Element>();
		for (Crack c : crackList) {
			Set<Element> splitElements = c.getSplitElements();
			Set<Element> crackTipElements = c.getCrackTipElements();

			for (Element e : splitElements) {
				for (Node n : e.getMeshPart().getNodes()) {
					if (enrichmentMap.get(n) == null) {
						enrichmentMap.put(n, new HashMap<Discontinuity, EnrichmentType>());
					}
					enrichmentMap.get(n).put(c, EnrichmentType.HEAVISIDE_ENRICHMENT);
				}
			}

			for (Element e : crackTipElements) {
				for (Node n : e.getMeshPart().getNodes()) {
					if (enrichmentMap.get(n) == null) {
						enrichmentMap.put(n, new HashMap<Discontinuity, EnrichmentType>());
					}
					enrichmentMap.get(n).put(c, EnrichmentType.CRACK_TIP_ENRICHMENT);
				}
			}

			fullyEnrichedElements.addAll(splitElements);
			fullyEnrichedElements.addAll(crackTipElements);
		}

		/*
		 * check for blending elements
		 */
		ArrayList<Element> blendingElements = new ArrayList<Element>();
		for (int i = 0; i < m.countElements(); i++) {

			Element e = m.getElement(i);

			if (!fullyEnrichedElements.contains(e)) {
				final MeshPart mp = e.getMeshPart();
				Node[] nodes = mp.getNodes();
				for (Node n : nodes) {
					if (enrichmentMap.get(n) != null) {
						blendingElements.add(e);
						break;
					}
				}
			}

		}

		/*
		 * enrich fully enrich elements
		 */
		for (Element e : fullyEnrichedElements) {
			final MeshPart mp = e.getMeshPart();
			Node[] nodes = mp.getNodes();

			FunctionBasisOnRnToR basis = BasisManager.getInstance().getBasis(mp.getType());

			FunctionBasisOnRnToR enrichedBasis = new FunctionBasisOnRnToR() {

				@Override
				public Object clone() {
					throw new Error("cloning is not supported for this class");
				}
			};

			/*
			 * check each node and define enrichment basis for that node
			 */
			for (int i = 0; i < nodes.length; i++) {
				final Node n = nodes[i];

				FunctionRnToR Ni = basis.getBasisFunction(i);
				enrichedBasis.addBasisFunction(Ni);

				int nModes = 1;
				Map<Discontinuity, EnrichmentType> enrichment = enrichmentMap.get(n);

				if (enrichment != null) {

					/*
					 * scan all cracks for one contains this enrichment
					 */
					for (final Crack c : crackList) {
						EnrichmentType enrichmentType = enrichment.get(c);
						final double[] xi = IntegrationManager.getInstance()
								.getReferenceCoordinate(mp.getType(), i);

						switch (enrichmentType) {

						case HEAVISIDE_ENRICHMENT:

							FunctionRnToR H = new FunctionRnToR() {

								@Override
								public double valueAt(double... x) {
									double psi = c.getLevelSet("psi").valueAt(mp, x);
									double psii = c.getLevelSet("psi").valueAt(n);
									return Math.signum(psi) - Math.signum(psii);
								}

								@Override
								public double[] gradientAt(double... x) {
									return new double[x.length];
								}
							};

							FunctionRnToR NH = new MultFunctionRnToR(1.0, Ni, H);
							enrichedBasis.addBasisFunction(NH);

							nModes += 1;

							break;

						case CRACK_TIP_ENRICHMENT:

							final FunctionRnToR[] B = BasisManager.getInstance()
									.createCrackTipBasis(mp, c);

							for (int j = 0; j < B.length; j++) {

								final int k = j;
								final double shift = B[k].valueAt(xi);
								// Debugger.watch("shift = ", shift);

								FunctionRnToR B_shift = new FunctionRnToR() {

									@Override
									public double valueAt(double... x) {
										return B[k].valueAt(x) - shift;
									}

									@Override
									public double[] gradientAt(double... x) {
										return B[k].gradientAt(x);
									}
								};

								enrichedBasis.addBasisFunction(new MultFunctionRnToR(1.0, Ni,
										B_shift));

							}

							nModes += 4;

							break;

						default:
							break;
						}
					}

					/*
					 * TODO: scan for other type of discontinuities
					 */

				}

				n.setNumberOfModes(nModes);

			}

			mp.setBasis(enrichedBasis);
		}

		/*
		 * enrich blending elements. Currently ramp function is not supported
		 */
		for (Element e : blendingElements) {
			final MeshPart mp = e.getMeshPart();
			Node[] nodes = mp.getNodes();

			FunctionBasisOnRnToR basis = BasisManager.getInstance().getBasis(mp.getType());

			FunctionBasisOnRnToR enrichedBasis = new FunctionBasisOnRnToR() {

				@Override
				public Object clone() {
					throw new Error("cloning is not supported for this class");
				}
			};

			/*
			 * add modified enrichment function
			 */
			for (int i = 0; i < nodes.length; i++) {
				final Node n = nodes[i];
				FunctionRnToR Ni = basis.getBasisFunction(i);

				enrichedBasis.addBasisFunction(Ni);

				Map<Discontinuity, EnrichmentType> enrichment = enrichmentMap.get(n);

				if (enrichment != null) {
					/*
					 * scan all cracks for one contains this enrichment
					 */
					for (int j = 0; j < crackList.size(); j++) {
						final Crack c = crackList.get(j);
						EnrichmentType enrichmentType = enrichment.get(c);
						final double[] xi = IntegrationManager.getInstance()
								.getReferenceCoordinate(mp.getType(), i);

						switch (enrichmentType) {

						case HEAVISIDE_ENRICHMENT:

							FunctionRnToR H = new FunctionRnToR() {

								@Override
								public double valueAt(double... x) {
									double psi = c.getLevelSet("psi").valueAt(mp, x);
									double psii = c.getLevelSet("psi").valueAt(n);
									return Math.signum(psi) - Math.signum(psii);
								}

								@Override
								public double[] gradientAt(double... x) {
									return new double[x.length];
								}
							};

							FunctionRnToR NH = new MultFunctionRnToR(1.0, Ni, H);
							enrichedBasis.addBasisFunction(NH);

							break;

						case CRACK_TIP_ENRICHMENT:

							final FunctionRnToR[] B = BasisManager.getInstance()
									.createCrackTipBasis(mp, c);

							for (int k = 0; k < B.length; k++) {

								final int l = k;
								final double shift = B[l].valueAt(xi);

								FunctionRnToR B_shift = new FunctionRnToR() {

									@Override
									public double valueAt(double... x) {
										return B[l].valueAt(x) - shift;
									}

									@Override
									public double[] gradientAt(double... x) {
										return B[l].gradientAt(x);
									}
								};

								MultFunctionRnToR NB = new MultFunctionRnToR(1.0, Ni, B_shift);
								enrichedBasis.addBasisFunction(NB);
							}

							break;

						default:
							break;
						}
					}
					/*
					 * TODO: scan for other type of discontinuities
					 */
				}
			}

			mp.setBasis(enrichedBasis);
		}

		/*
		 * for non-enrich element, switch back to normal isoparametric basis
		 */
		for (int i = 0; i < m.countElements(); i++) {
			Element e = m.getElement(i);
			if (!fullyEnrichedElements.contains(e) && !blendingElements.contains(e)) {
				MeshPart mp = e.getMeshPart();
				mp.setBasis(BasisManager.getInstance().getBasis(mp.getType()));
			}
		}

		/***********************************************************************
		 * build integration rule for enriched elements. Currently, this code
		 * does not support for crack intersection
		 ***********************************************************************/
		IntegrationManager im = IntegrationManager.getInstance();

		/*
		 * switch all elements to normal integration rule
		 */
		for (int i = 0; i < m.countElements(); i++) {
			Element e = m.getElement(i);
			e.setIntegrationOrder(im.getDefaultIntegrationOrder(e.getMeshPart().getType()));
		}

		for (Crack c : crackList) {

			/**********************************
			 * processing the split elements
			 **********************************/
			Set<Element> splitElements = c.getSplitElements();

			for (Element e : splitElements) {
				MeshPart mp = e.getMeshPart();
				int dim = mp.getDimension();
				Node[] nodes = mp.getNodes();

				/*
				 * separate nodes into two domains
				 */
				ArrayList<double[]> posDomain = new ArrayList<double[]>();
				ArrayList<double[]> negDomain = new ArrayList<double[]>();
				for (int i = 0; i < nodes.length; i++) {
					double[] xi = im.getReferenceCoordinate(mp.getType(), i);
					if (c.getLevelSet("psi").valueAt(nodes[i]) > 0.0) {
						posDomain.add(xi);
					} else {
						negDomain.add(xi);
					}
				}

				/*
				 * compute intersection points
				 */
				ArrayList<double[]> intersects = c.computeIntersectionAt(e);
				// Debugger.watch(e.getId());
				// Debugger.watch("", intersects);

				if (intersects.size() != 2) {
					Debugger.watch("intersection: ", intersects);
					Debugger.watch("element id = " + e.getId());
					throw new Error("number of intersection points is not 2");
				}

				/*
				 * detect the kinks/high curvature points. TODO
				 */

				/*
				 * generate center points in each domain
				 */
				double[] p1 = intersects.get(0).clone();
				LinAlg.add(dim, 1.0, intersects.get(1), 1.0, p1);
				int cnt = 2;
				for (double[] p : posDomain) {
					LinAlg.add(dim, 1.0, p, 1.0, p1);
					cnt++;
				}
				LinAlg.scale(dim, 1.0 / cnt, p1);
				posDomain.add(p1);

				double[] p2 = intersects.get(0).clone();
				LinAlg.add(dim, 1.0, intersects.get(1), 1.0, p2);
				cnt = 2;
				for (double[] p : negDomain) {
					LinAlg.add(dim, 1.0, p, 1.0, p2);
					cnt++;
				}
				LinAlg.scale(dim, 1.0 / cnt, p2);
				negDomain.add(p2);

				/*
				 * build the Delaunay triangulation
				 */
				AbstractDelaunay d1 = new Delaunay2(-2, 2, -2, 2);
				AbstractDelaunay d2 = new Delaunay2(-2, 2, -2, 2);
				for (double[] p : posDomain) {
					d1.addPoint(p[0], p[1]);
				}
				for (double[] p : negDomain) {
					d2.addPoint(p[0], p[1]);
				}
				List<Triangle> triangles = d1.getTriangles();
				triangles.addAll(d2.getTriangles());

				/*
				 * generate integration points for each domain
				 */
				ArrayList<IntegrationPoint> ips = im.createIntegrationPoints(p_on_sub_triangle,
						triangles);
				im.setIntegrationPoints(e, ips);

				/*
				 * change the integration rule
				 */
				e.setIntegrationType(IntegrationType.HEAVISIDE_SPLIT_INTEGRATION);

				/*
				 * TODO:mapping variables at existing integration points to new
				 * integration points
				 */
				e.Initialize();
			}

			/**********************************
			 * processing the crack tip elements
			 **********************************/
			Set<Element> crackTipElements = c.getCrackTipElements();
			for (Element e : crackTipElements) {
				MeshPart mp = e.getMeshPart();
				Node[] nodes = mp.getNodes();

				/*
				 * compute intersection points
				 */
				ArrayList<double[]> intersects = c.computeIntersectionAt(e);
				assert intersects.size() == 2 : "something wrong with the intersection";

				/*
				 * compute crack tip in reference coordinate
				 */
				double[] crackTip = c.computeLevelSetIntersection(e);

				/*
				 * build the Delaunay triangulation
				 */
				AbstractDelaunay d = new Delaunay2(-2, 2, -2, 2);
				for (int i = 0; i < nodes.length; i++) {
					double[] xi = im.getReferenceCoordinate(mp.getType(), i);
					d.addPoint(xi[0], xi[1]);
				}
				for (double[] p : intersects) {
					d.addPoint(p[0], p[1]);
				}
				d.addPoint(crackTip[0], crackTip[1]);

				// Viewer v = new Viewer();
				// d.plot(v);
				// v.setVisible(true);

				/*
				 * generate integration points
				 */
				ArrayList<IntegrationPoint> ips = im.createIntegrationPoints(p_on_sub_triangle,
						d.getTriangles());
				im.setIntegrationPoints(e, ips);

				/*
				 * change the integration rule
				 */
				e.setIntegrationType(IntegrationType.CRACK_TIP_INTEGRATION);

				/*
				 * TODO:mapping variables at existing integration points to new
				 * integration points
				 */
				e.Initialize();

			}

			Debugger.watch("split elements");
			for (Element e : splitElements) {
				Debugger.watch(e.getId());
			}

			Debugger.watch("crack tip elements");
			for (Element e : crackTipElements) {
				Debugger.watch(e.getId());
			}
		}

		/*
		 * set the integration order for blending elements TODO: check the order
		 */
		for (Element e : blendingElements) {
			e.setIntegrationOrder(2 * e.getIntegrationOrder());
			e.Initialize();
		}

		/*
		 * re-enumerate the mesh
		 */
		m.preAnalysis();

		Debugger.watch("number of elements: " + m.countElements());
		Debugger.watch("enriched elements: ");
		for (Element e : fullyEnrichedElements) {
			System.out.print(e.getId() + ", ");
		}
		System.out.println();
		Debugger.watch("blending elements: ");
		for (Element e : blendingElements) {
			System.out.print(e.getId() + ", ");
		}
		System.out.println();

	}

	@Override
	public void fireIterationStarted() {
	}

	@Override
	public void fireIterationFinished() {
	}

	@Override
	public void fireTimeStepFinished() {
		for (Crack c : crackList) {
			/*
			 * compute stress intensity factor
			 */
			double[] sifs = c.computeSIFs(m);
			double K1 = sifs[0];
			double K2 = sifs[1];
			double thetac = 2 * Math.atan(0.25 * (K1 / K2 - Math.signum(K2)
					* Math.sqrt(Math.pow(K1 / K2, 2) + 8)));
			/*
			 * update level set
			 */
			c.propagate(c.getPropagationStep(), thetac);
		}
	}

	@Override
	public void fireFinished() {

	}

}
